DSA 2021 - TD R
Fabien FAIVRE
18 juillet 2021
1 Présentation du projet
Ce projet s’inscrit dans le cadre de la formation DSA, et évalue l’apprentissage de R réalisé auprès de Robin RYDER.
Le thème de ce projet est libre.
J’ai choisi d’étudier la manière dont les approches de machine learning (ci après
ML) peuvent compléter les approches classiques par modèles linéaires généralisés (GLM) dans le cadre de la tarification en assurance non-vie.
1.1 Les bonnes propriétés du GLM dans le cadre d’une modélisation de la fréquence de sinistres
Soit \(N_{i}\) le nombre de sinistre de la police \({i}\), alors dans le modèle Poisson homogène on suppose \(N_{i} \sim \mathscr{P}(\lambda \nu_{i})\) où \(\nu_{i}\) est l’exposition au risque de la police \(i\), \(\lambda\) est la fréquence de sinistre annuelle (homogène sur le portefeuille)et \(\mathscr{P}\) la loi de Poisson. On pose \(Y_{i} = \frac{N_{i}}{\nu_{i}}\) la fréquence annuelle des sinistres pour la police \(i\) et \(\mathbf{x_{i}}=(x_{1},x_{2},..., x_{q})\) les \(q\) caractéristiques de la police \(i\).
En tarification, on suppose dans les faits que \(\lambda\) dépend de \(\mathbf{x}\), de sorte que \(\mathbb{P}[Y=k/\nu] = \mathbb{P}[N=k] = e^{-\lambda(\mathbf{x})\nu} \frac{(\lambda(\mathbf{x})\nu)^{k}}{k!}\) avec \(ln(\lambda(\mathbf{x})) = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + \cdots + \beta_{q}x_{q} \overset{def}{=} \langle \mathbf{\beta}, \mathbf{x} \rangle\) en complétant \(\mathbf{x}\) par \(x_{0} = 1\).
Lorsqu’on utilise les approches classiques de GLM, on estime \(\hat{\lambda}\) en utilisant l’estimateur du maximum de vraisemblance \(\hat{\beta}\) de \(\beta\).
Cet estimateur peut être obtenu soit en maximisant la log vraisemblance, soit de manière équivalente en minimisant la déviance de Poisson.
En complément, lorsqu’on utilise un GLM incluant un intercept et sa fonction de lien canonique, l’utilisatiuon d’un estimateur du maximum de vraisemblance apporte de bonnes propriétés pour son usage en tarification dont celle dite de propriété d’équilibre (Property 2.4 p33 de Data Analytics for Non-Life Insurance Pricing), qui précise que :
\(\displaystyle\sum_{i=1}^{n} \nu_i\hat{\lambda}(\mathbf{x_{i}}) = \displaystyle\sum_{i=1}^{n} \nu_i exp\langle \mathbf{\hat{\beta}}, \mathbf{x_{i}} \rangle = \displaystyle\sum_{i=1}^{n} N_{i}\)
De sorte que le nombre total de sinistres estimés est égale au nombre total de sinistres observés.
Cette propriété est même plus étendue dans certains modèles GLM, notamment lorsque le modèle poissonien avec lien log (exemple (b) p423 de A systematic relationship between minimum bias and gener- alized linear models) : dans ce cas l’équilibre existe au global du portefeuille, mais aussi pour chaque classe des variables catégorielles incluses dans le modèle.
Ces propriétés d’équilibres sont critiques pour un assureur dans un exercice de tarification, puisqu’elles assurent qu’une structure tarifaire basée sur cette approche lui permettra d’équilibrer (au moins théoriquement) ses encaissement avec les sinistres qu’il s’est engagé à indemniser.
Le GLM est donc interprétable, souple à l’usage et présente des propriétés d’équilibre global et local. Néanmoins son ajustement nécessite un temps humain important et le praticien doit s’en remettre à son expérience pour estimer les éventuelles interactions à inclure dans son analyse.
1.2 L’apport et les limites des modèles ML en tarification
Les modèles ML par leur grande flexibilité sont des candidats naturels pour remplacer ou à minima compléter les modèles GLM du fait de leur performance sur une grande variété de cas d’usages ces dernières années. Toutefois, ces modèles nécessitent qulques ajustements avant de pouvoir être utilisés dans des problèmes de tarification.
1.2.1 Ajustement #1 : définition de la fonction de perte
Une première adaptation des modèles classiques consiste donc à définir la déviance de poisson comme fonction de perte des algorithmes de ML.
C’est notamment l’approche préconisée par dans Data Analytics for Non-Life Insurance Pricing de Mario V. Wuthrich et Christoph Buser.
1.2.2 Limites de la modification de la fonction de perte
La définition d’une fonction de perte appropriée est devenue l’approche par défaut de l’utilisation des méthodes ML en tarification. Néanmoins, dès 2019 Wütrich remarqua que les sinistres estimés par les modèles ML ne correspondaient pas à la sinistralité observée sur le portefeuille, même sur le jeu d’entraînement
Il semble que l’utilisation directe de modèles de ML, mêmes adaptés pour minimiser la déviance, ne présentent pas les garanties d’équilibre des GLM.
Les modèles de ML permettent une meilleurs corrélation avec la réponse individuelle de chaque police, du fait de la plus grande flexibilité accordée à la forme de la contribution (éventuellement non linéaire) de chaque variable à la sortie. Toutefois, cette meilleure corrélation s’obtient au prix de l’abandon de la notion d’équilibre. En effet, aucune garantie n’est expressément introduite au global du portefeuille ce qui peut produire des divergences, potentiellement importantes, entre la somme des primes demandées et les sinistres que l’assureur s’est engagé à régler.
Cet effet est attribué par M Wütrich aux conditions d’arrêt précoce des algorithmes d’optimisation utilisés par les modèles ML.
1.3 Pistes de réconciliation des approches GLM et ML pour la tarification
Deux options s’offrent dès lors à l’actuaire :
essayer d’extraire des informations de la manière dont les modèles ML ont appris les efefts des variables sur la sortie pour compléter les GLM classiquement manipulés avec l’objectiof de combler en partie l’écart de performance (s’il est important) avec le modèle ML tout en conservant les propriiétés et l’interprétabilité naturelle des modèles GLM. C’est la voie suivie par l’article Peeking into the black box de Christian Lorentzen et Michael Mayer
corriger les approches ML en les combinant avec une approche GLM pour en récupérér les propriétés. Le gain de cette approche devrait être plus important que pour la méthode précédente tout en présentant les propriétés d’équilibre des modèles ML sans toutefois apporter l’interprétabilité de ces modèles. C’est la proposition de de Michel Denuit, Arthur Charpentier et Julien Trufin dans Autocalibration and Tweedie-dominance for Insurance Pricing with Machine Learning repartent du constat de Wütrich et proposent une solution de contournement sous la forme d’un ajustement en deux temps successifs :
Dans un premier temps, un modèle ML est ajusté pour déterminer une première estimation \(\hat{\pi}\) des primes théoriques
Dans un deuxième temps, ces primes théoriques \(\hat{\pi}\) sont utilisées comme variable entrante d’un modèle d’autocalibration local pour produire des primes corrigées \(\widehat{\pi_{BC}}\). Cette méthode consiste à ajuster des modèles GLM locaux réduits à un intercept sur des ensembles de points “proches”. La notion de proximité ici se fait qui sont définis à partir de la proximité des primes théoriques estimées par ML \(\hat{\pi}\) (l’étape de ML rassemble les observations qui se ressemblent, l’étape locale GLM lisse les écarts entre observations voisines pour retrouver les propriétés d’équilibre).
1.4 Objectif du projet
Dans ce rapport, faute de temps, nous n’explorerons que l’option 1 en commençant par comparer la performance d’un modèle GLM classique avec un modèle GBM et un modèle XGBoost.
Nous tâcherons de comprendre la manière dont les modèles ML ont appris les relations entre les variables pour dans un troisième temps proposer une version améliorée du GLM initial.
Pour ce faire, nous prendrons comme base d’analyse le jeu de données tarifaire d’un portefeuille MRH qui nous avait été confié lors du Hackathon.
L’enjeu étant surtout d’illustrer le principe, nous nous restreindrons dans la suite de ce document à l’analyse de la seule fréquence des sinistres.
2 Analyse
2.1 Mise en place de l’environnement
Le code suivant est un ensemble d’utilitaire pour naviguer dans les répertoires relatifs du projet et installer ou charger les packages nécessaires à l’analyse :
if (!require("here")){
install.packages("here")
library("here")
}
set_here
LoadPackage <- function (load.lib=c("")) {
install.lib<-load.lib[!load.lib %in% installed.packages()]
for(lib in install.lib) install.packages(lib,dependencies=TRUE)
sapply(load.lib,require,character=TRUE)
}
LoadPackage(c('ggplot2', 'dplyr', 'formattable', 'DT', 'tidyr', 'caret', 'plotly', 'xgboost', 'flashlight', 'MetricsWeighted', 'corrplot', 'lsr', 'h2o', 'wesanderson', 'locfit'))2.2 code pour la phase exploratoire
Dans cette section nous allons définr deux fonctions qui seront utilisées de manière répétées dans la phase exploratoire du jeu de données à savoir une fonction générique de tracé des exprositions et fréquences moyennes de sinistres d’une part et de création de tables synthétiques d’autre part :
plot_obs<-function(df,feature,use_year=T) {
if (use_year) {
dt <-df %>%
group_by(df[[feature]], ANNEE, .drop=FALSE) %>%
summarise(Freq = sum(NB*EXPO)/sum(EXPO), nb_lignes = n(), EXPO=sum(EXPO) , pct_EXPO = sum(EXPO) / sum(df$EXPO)) %>%
rename(feat = 'df[[feature]]')
} else {
dt <-df %>%
group_by(df[[feature]], .drop=FALSE) %>%
summarise(Freq = sum(NB*EXPO)/sum(EXPO), nb_lignes = n(), EXPO=sum(EXPO) , pct_EXPO = sum(EXPO) / sum(df$EXPO)) %>%
rename(feat = 'df[[feature]]')
}
fig <-plot_ly()
fig <- fig %>% layout(title = paste0('Fréquence de sinistres par ', feature),
xaxis = list(title = feature),
yaxis = list(side = 'left', title = 'Exposition', showgrid = FALSE, zeroline = FALSE),
yaxis2 = list(side = 'right', overlaying = "y", title = 'Fréquence', showgrid = FALSE, zeroline = FALSE))
if (use_year) {
for (year in unique(dt$ANNEE)) {
dt_temp <- dt %>% filter(ANNEE==year)
fig <- fig %>% add_trace(
x= dt_temp$feat, y= dt_temp$EXPO, type = 'bar', name = paste0('Exposition_', year),
hovertemplate = 'Expo: %{y:.0f}<br>'
)
fig<- fig %>% add_trace(
x= dt_temp$feat, y= dt_temp$Freq, type = 'scatter', mode = 'lines', yaxis = 'y2', name= paste0('Fréquence_', year),
hovertemplate = 'Freq: %{y:.2%}<br>'
)
}
} else {
dt_temp <- dt
fig <- fig %>% add_trace(
x= dt_temp$feat, y= dt_temp$EXPO, type = 'bar', name = 'Exposition',
hovertemplate = 'Expo: %{y:.0f}<br>'
)
fig<- fig %>% add_trace(
x= dt_temp$feat, y= dt_temp$Freq, type = 'scatter', mode = 'lines', yaxis = 'y2', name= 'Fréquence',
hovertemplate = 'Freq: %{y:.2%}<br>'
)
}
return(fig)
}resume <- function(df, feature) {
dt <-df %>%
group_by(df[[feature]], .drop=FALSE) %>%
summarise(Freq = sum(NB*EXPO)/sum(EXPO), nb_lignes = n(), EXPO=sum(EXPO) , pct_EXPO = sum(EXPO) / sum(df$EXPO ))
colnames(dt)[1] <- feature
table<- formattable(dt,
align = c("l", rep("r", NCOL(dt) - 1)),
list(
Freq = percent,
nb_lignes= accounting,
EXPO = accounting,
pct_EXPO = percent
)
)
fig1<- plot_obs(df, feature, use_year = F)
fig2<- plot_obs(df, feature, use_year = T)
return(list(dt=dt, table=table, fig=fig1, figyr=fig2))
}3 Présentation du jeu de données
Le jeu de données est constitué de 4 fichiers :
list.files(path=here('data', 'raw'))## [1] "Data" "Data.zip" "dsa2021_hackathon.zip"
## [4] "expo_test.csv" "expo_train.csv" "primes2019.csv"
## [7] "sin_train.csv"Le projet suit un portefeuille d’assurance MRH suivi sur plusieurs années de 2016 à 2018 inclus. L’objectif initial du Hackathon est d’estimer les primes 2019 sur un échantillon de contrats.
La nature de la garantie modélisée ne nous a pas été communiquée.
Pour éviter les redondances, les analyses descriptives sont concentrées sur le fichier combiné en fin de section.
3.1 expo_train
Le fichier expo_train.csv contient la liste des contrats en risques historiquement suivis :
expo_train = read.table(file = here('data', 'raw', 'expo_train.csv'), header=T, sep=',', dec='.', encoding = 'UTF-8', stringsAsFactors = T)
datatable(head(expo_train))str(expo_train)## 'data.frame': 155651 obs. of 19 variables:
## $ X : int 384538 441079 119668 5124 130065 450830 151401 296526 71801 360445 ...
## $ EXPO : num 1 0.825 1 1 1 ...
## $ FORMULE : Factor w/ 4 levels "ALL_INCLUDE",..: 4 2 3 3 3 2 3 2 4 2 ...
## $ TYPE_RESIDENCE : Factor w/ 2 levels "PRINCIPALE","SECONDAIRE": 1 1 1 2 1 1 1 1 1 2 ...
## $ TYPE_HABITATION : Factor w/ 2 levels "APPARTEMENT",..: 1 2 1 2 2 2 2 2 2 1 ...
## $ NB_PIECES : int 1 NA 3 2 1 2 2 2 1 3 ...
## $ SITUATION_JURIDIQUE: Factor w/ 2 levels "LOCATAIRE","PROPRIO": 2 2 1 1 1 1 1 1 2 1 ...
## $ NIVEAU_JURIDIQUE : Factor w/ 2 levels "JUR1","JUR2": 1 1 1 1 1 1 1 1 1 1 ...
## $ VALEUR_DES_BIENS : num 3500 0 35000 9000 20000 9000 20000 9000 3500 0 ...
## $ OBJETS_DE_VALEUR : Factor w/ 2 levels "NIVEAU_1","NIVEAU_2": 1 1 1 1 1 1 1 1 1 1 ...
## $ ZONIER : Factor w/ 82 levels "A1","A10","A11",..: 48 3 39 75 82 71 70 81 14 12 ...
## $ NBSIN_TYPE1_AN1 : int 0 0 0 0 0 0 1 0 0 0 ...
## $ NBSIN_TYPE1_AN2 : int 1 1 1 1 1 1 1 1 1 1 ...
## $ NBSIN_TYPE1_AN3 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ NBSIN_TYPE2_AN1 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ NBSIN_TYPE2_AN2 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ NBSIN_TYPE2_AN3 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ id : int 5 9 11 13 14 15 20 21 28 30 ...
## $ ANNEE : int 2017 2018 2017 2018 2017 2017 2018 2016 2018 2016 ...Les données ont été prétraitées pour disposer d’une ligne par période de risque annuelle homogène.
Pour des raisons d’anonymisation, les identifiants contrats ont été supprimés de sorte qu’il n’est pas possible de suivre l’évolution du risque d’une période sur l’autre. Pour la même raison, il ne sera pas possible dans la construction des jeux de validation et de test de s’assurer que toute l’expérience d’un même assuré est bien intégralement dans un jeu d’entraînement, de validation ou de test.
La signification des colonnes présentes est la suivante :
X: variable non nommée. Il s’agit d’un artefact du processus de création du fichier initial. On ne la prend pas en compteEXPO: exposition en année risque du contrat. Sa valeur précalculée pour chaque ligne est comprise entre 0 et 1FORMULE: variable catégorielle codant la formule du contrat comprenant 3 niveauxMEDIUM,ESSENTIELetCONFORTTYPE_RESIDENCE: variable catégorielle codant le fait que le bien est une résidencePRINCIPALE, ouSECONDAIRETYPE_HABITATION: variable catégorielle codant le fait que le bien est unAPPARTEMENT, ou uneMAISONNB_PIECES: variable numérique indiquant le nombre de pièces du logementSITUATION_JURIDIQUE: variable catégorielle indiquant si le souscripteur est prorpiétaire (PROPRIO) ou locataire (LOCATAIRE) du logement assuréNIVEAU_JURIDIQUE: variable catégorielle codant le niveau de couverture de la garantie juridiqueVALEUR_DES_BIENS: variable numérique reflétant la valeur couverte du contenu du logementOBJETS_DE_VALEUR: variable catégorielle codant le niveau de couverture des objetsZONIER: variable catégorielle codant la zone du bien assuré. Le code est constitué d’une lettre représentant une zone et d’un nombre représentant une partie de la zoneNBSIN_TYPE1_AN1: variable numérique indiquant le nombre de sinistre de type 1 l’année précédenteNBSIN_TYPE1_AN2: variable numérique indiquant le nombre de sinistre de type 1 il y a 2 ansNBSIN_TYPE1_AN3: variable numérique indiquant le nombre de sinistre de type 1 il y a 3 ansNBSIN_TYPE2_AN1: variable numérique indiquant le nombre de sinistre de type 2 l’année précédenteNBSIN_TYPE2_AN2: variable numérique indiquant le nombre de sinistre de type 2 il y a 2 ansNBSIN_TYPE2_AN3: variable numérique indiquant le nombre de sinistre de type 2 il y a 3 ansid: identifiant du risque, clef de jointureANNEE: variable catégorielle indiquant l’année d’observation du risque :2016,2017,2018
3.2 sin_train
Le fichier expo_train.csv contient la liste des contrats en risques historiquement suivis :
sin_train = read.table(file = here('data', 'raw', 'sin_train.csv'), header=T, sep=';', dec=',', encoding = 'UTF-8', stringsAsFactors = T)
datatable(head(sin_train))str(sin_train)## 'data.frame': 2693 obs. of 4 variables:
## $ id : int 15 277 643 668 730 1816 1852 2061 2270 2322 ...
## $ NB : int 1 1 1 1 1 1 1 1 1 1 ...
## $ COUT : num 521.4 3000.2 26.3 462.4 640.9 ...
## $ ANNEE: int 2017 2016 2018 2017 2018 2018 2016 2018 2018 2017 ...id: identifiant du risque, clef de jointure avec le fichierexpo_trainNB: Nombre de sinistres de la typologie modélisée (inconnue) sur la période d’observation
| NB | nb_lignes |
|---|---|
| 1 | 2,613.00 |
| 2 | 79.00 |
| 3 | 1.00 |
3.3 fichier combiné
On peut désormais combine les information d’exposition et de sinistres :
mrh <- expo_train %>%
left_join(sin_train, by =c('id','ANNEE')) %>%
replace_na(list('NB'=0, 'COUT'=0))N’ayant pas accès au fichier de test utilisé lors du Hackathon, on dissocie ici le fichier mrh en un fichier train (70%), un fichier val (20%) et un fichier test (10%). Cette approche est prise ici sans considération des années d’étude.
C’est un biais fort qui aurait été rédhibitoire dans la réalité puisqu’il aurait été de nature à induire des effets de leakage entre les jeu d’entrainement et de test.
Toutefois cette approche a été prise par commodité et pour garder une proportion significative de volume dans les jeux de validation et test et parce que l’enjeu de l’analyse n’est pas nécessairement la performance réelle des algorithmes, mais plutôt la démarche pour augmenter les modèles GLM à partir des informations extraites des modèles ML.
Nous utilisons ici la fonction createDataPartition du package caret qui permet de faire du stratified sampling par rapport à une variable d’entrée.
set.seed(42)
trainIndex<- createDataPartition(mrh$NB>=0, p=.7, list=FALSE, time=1)
mrh_train<-mrh[trainIndex,]
mrh_test<-mrh[ -trainIndex,]
valIndex<- createDataPartition(mrh_test$NB>=0, p=.66, list=FALSE, time=1)
mrh_val<-mrh_test[ valIndex,]
mrh_test<-mrh_test[ -valIndex,]3.3.1 Taux de sinistres moyen :
tx <- sum(mrh_train$NB)/sum(mrh_train$EXPO)
print(paste0('taux de sinistre moyen annuel : ', round(tx * 100,2), '%'))## [1] "taux de sinistre moyen annuel : 2.12%"3.3.2 Taux de sinistre par variable :
3.3.2.1 FORMULE
| FORMULE | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| ALL_INCLUDE | 0.57% | 4,695.00 | 3,355.67 | 3.67% |
| CONFORT | 2.11% | 61,372.00 | 51,858.85 | 56.79% |
| ESSENTIEL | 2.22% | 15,555.00 | 13,597.27 | 14.89% |
| MEDIUM | 1.69% | 27,334.00 | 22,501.75 | 24.64% |
ENSEIGNEMENT : Le comportement observé est stable dans le temps. Les niveaux CONFORT et ESSENTIEL semblent avoir un niveau de risque de fréquence équivalent
3.3.2.2 TYPE_RESIDENCE
| TYPE_RESIDENCE | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| PRINCIPALE | 2.15% | 91,410.00 | 78,937.92 | 86.45% |
| SECONDAIRE | 0.75% | 17,546.00 | 12,375.63 | 13.55% |
ENSEIGNEMENT : Le comportement observé est stable dans le temps. Les résidences SECONDAIRE ont systématiquement une fréquence moindre
3.3.2.3 TYPE_HABITATION
| TYPE_HABITATION | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| APPARTEMENT | 1.98% | 49,749.00 | 43,093.59 | 47.19% |
| MAISON | 1.95% | 59,207.00 | 48,219.96 | 52.81% |
ENSEIGNEMENT : Le comportement observé N’est PAS stable dans le temps. Le type d’habitation est à exclure.
3.3.2.4 NB_PIECES
| NB_PIECES | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| 0 | 0.12% | 2,610.00 | 1,715.52 | 1.88% |
| 1 | 1.61% | 28,256.00 | 22,425.90 | 24.56% |
| 2 | 1.83% | 42,724.00 | 36,390.78 | 39.85% |
| 3 | 2.48% | 22,181.00 | 19,440.34 | 21.29% |
| 4 | 2.95% | 6,639.00 | 5,879.31 | 6.44% |
| NA | 1.98% | 6,546.00 | 5,461.69 | 5.98% |
ENSEIGNEMENT : Le NB_PIECES présente des valeurs manquantes. La fréquence tend à augmenter avec le nombre de pièces.
3.3.2.5 SITUATION_JURIDIQUE
| SITUATION_JURIDIQUE | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| LOCATAIRE | 2.21% | 60,339.00 | 53,225.98 | 58.29% |
| PROPRIO | 1.63% | 48,617.00 | 38,087.57 | 41.71% |
ENSEIGNEMENT : Les propriétaires semblent avoir un fréquence moindre que les locataires de manière consistante, mais avec de la variabilité.
3.3.2.6 NIVEAU_JURIDIQUE
| NIVEAU_JURIDIQUE | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| JUR1 | 1.97% | 107,497.00 | 90,273.20 | 98.86% |
| JUR2 | 1.52% | 1,459.00 | 1,040.35 | 1.14% |
ENSEIGNEMENT : Le NIVEAU_JURIDIQUE N’est PAS stable dans le temps. A exclure.
3.3.2.7 VALEUR_DES_BIENS
| VALEUR_DES_BIENS | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| 0 | 0.25% | 5,782.00 | 3,682.96 | 4.03% |
| 3500 | 1.36% | 36,543.00 | 28,888.32 | 31.64% |
| 9000 | 1.95% | 24,283.00 | 21,014.51 | 23.01% |
| 20000 | 2.25% | 23,639.00 | 20,768.83 | 22.74% |
| 35000 | 2.86% | 8,047.00 | 7,187.60 | 7.87% |
| 50000 | 3.14% | 5,470.00 | 4,991.12 | 5.47% |
| 80000 | 3.09% | 4,557.00 | 4,210.03 | 4.61% |
| 100000 | 3.72% | 635.00 | 570.18 | 0.62% |
ENSEIGNEMENT : La fréquence augmente globalement avec la valeur des biens selon une relation non linéaire Par ailleurs au delà de 50k, l’expérience n’est pas stable. Les modalités supérieures à 80K devraient être groupées.
3.3.2.8 OBJETS_DE_VALEUR
| OBJETS_DE_VALEUR | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| NIVEAU_1 | 1.84% | 102,196.00 | 85,308.12 | 93.42% |
| NIVEAU_2 | 3.77% | 6,760.00 | 6,005.44 | 6.58% |
ENSEIGNEMENT : Le NIVEAU_2 présente une fréquence significativement plus élevée.
3.3.2.9 ZONIER
| ZONIER | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| A1 | 1.75% | 1,271.00 | 1,098.88 | 1.20% |
| A10 | 1.21% | 768.00 | 659.36 | 0.72% |
| A11 | 0.62% | 3,349.00 | 2,901.64 | 3.18% |
| A12 | 1.10% | 2,182.00 | 1,896.53 | 2.08% |
| A13 | 0.97% | 1,326.00 | 1,112.93 | 1.22% |
| A14 | 1.81% | 596.00 | 497.59 | 0.54% |
| A2 | 1.14% | 1,199.00 | 1,004.07 | 1.10% |
| A3 | 0.57% | 814.00 | 705.81 | 0.77% |
| A4 | 0.00% | 55.00 | 40.92 | 0.04% |
| A5 | 0.00% | 11.00 | 10.24 | 0.01% |
| A6 | 0.72% | 495.00 | 427.35 | 0.47% |
| A7 | 0.69% | 170.00 | 145.53 | 0.16% |
| A8 | 2.29% | 48.00 | 43.73 | 0.05% |
| A9 | 0.56% | 2,885.00 | 2,488.88 | 2.73% |
| B1 | 0.41% | 556.00 | 492.90 | 0.54% |
| B10 | 1.72% | 2,492.00 | 2,140.03 | 2.34% |
| B11 | 0.00% | 6.00 | 4.80 | 0.01% |
| B12 | 2.37% | 249.00 | 210.74 | 0.23% |
| B13 | 6.23% | 37.00 | 32.09 | 0.04% |
| B14 | 0.00% | 83.00 | 69.97 | 0.08% |
| B16 | 0.76% | 147.00 | 130.79 | 0.14% |
| B17 | 1.79% | 2,442.00 | 2,071.38 | 2.27% |
| B18 | 0.00% | 6.00 | 4.25 | 0.00% |
| B19 | 2.42% | 672.00 | 572.48 | 0.63% |
| B2 | 1.09% | 308.00 | 274.89 | 0.30% |
| B20 | 0.00% | 6.00 | 5.58 | 0.01% |
| B21 | 0.00% | 30.00 | 25.92 | 0.03% |
| B22 | 0.00% | 41.00 | 36.37 | 0.04% |
| B23 | 0.78% | 150.00 | 128.79 | 0.14% |
| B24 | 2.56% | 91.00 | 78.21 | 0.09% |
| B25 | 1.33% | 436.00 | 375.81 | 0.41% |
| B26 | 1.24% | 1,031.00 | 884.64 | 0.97% |
| B27 | 0.84% | 133.00 | 119.69 | 0.13% |
| B28 | 1.35% | 173.00 | 148.10 | 0.16% |
| B29 | 1.83% | 2,979.00 | 2,541.22 | 2.78% |
| B3 | 0.75% | 1,394.00 | 1,210.68 | 1.33% |
| B30 | 0.96% | 355.00 | 310.94 | 0.34% |
| B31 | 0.00% | 56.00 | 51.74 | 0.06% |
| B32 | 1.81% | 4,913.00 | 4,134.65 | 4.53% |
| B33 | 0.85% | 134.00 | 117.90 | 0.13% |
| B34 | 1.23% | 518.00 | 446.21 | 0.49% |
| B35 | 0.00% | 3.00 | 2.28 | 0.00% |
| B36 | 1.33% | 319.00 | 266.18 | 0.29% |
| B37 | 0.00% | 108.00 | 94.32 | 0.10% |
| B38 | 5.38% | 67.00 | 55.78 | 0.06% |
| B39 | 1.37% | 1,804.00 | 1,542.57 | 1.69% |
| B4 | 0.73% | 1,324.00 | 1,134.45 | 1.24% |
| B40 | 2.41% | 4,045.00 | 3,427.84 | 3.75% |
| B41 | 2.79% | 282.00 | 242.19 | 0.27% |
| B42 | 0.00% | 22.00 | 20.55 | 0.02% |
| B43 | 1.72% | 4,229.00 | 3,553.62 | 3.89% |
| B44 | 0.00% | 5.00 | 4.09 | 0.00% |
| B45 | 0.83% | 317.00 | 268.14 | 0.29% |
| B5 | 0.00% | 24.00 | 20.47 | 0.02% |
| B6 | 1.30% | 91.00 | 77.15 | 0.08% |
| B7 | 2.84% | 243.00 | 211.00 | 0.23% |
| B8 | 0.00% | 152.00 | 136.97 | 0.15% |
| B9 | 1.93% | 437.00 | 374.18 | 0.41% |
| C1 | 2.68% | 1,923.00 | 1,619.36 | 1.77% |
| C10 | 2.92% | 638.00 | 515.30 | 0.56% |
| C11 | 3.19% | 906.00 | 739.61 | 0.81% |
| C12 | 3.71% | 1,700.00 | 1,335.15 | 1.46% |
| C13 | 2.09% | 882.00 | 681.10 | 0.75% |
| C14 | 2.77% | 962.00 | 769.95 | 0.84% |
| C15 | 3.70% | 851.00 | 666.75 | 0.73% |
| C16 | 2.56% | 894.00 | 707.78 | 0.78% |
| C17 | 2.18% | 592.00 | 463.76 | 0.51% |
| C18 | 1.85% | 526.00 | 429.41 | 0.47% |
| C19 | 3.29% | 2,006.00 | 1,607.93 | 1.76% |
| C2 | 2.34% | 4,216.00 | 3,597.31 | 3.94% |
| C20 | 3.06% | 10,231.00 | 8,576.45 | 9.39% |
| C21 | 1.62% | 3,365.00 | 2,789.43 | 3.05% |
| C22 | 1.50% | 1,789.00 | 1,430.49 | 1.57% |
| C23 | 2.15% | 4,418.00 | 3,642.23 | 3.99% |
| C24 | 1.91% | 1,360.00 | 1,122.90 | 1.23% |
| C3 | 1.18% | 2,127.00 | 1,817.00 | 1.99% |
| C4 | 3.82% | 1,174.00 | 974.79 | 1.07% |
| C5 | 1.79% | 4,772.00 | 3,952.35 | 4.33% |
| C6 | 1.59% | 3,614.00 | 3,073.13 | 3.37% |
| C7 | 2.10% | 1,442.00 | 1,195.34 | 1.31% |
| C8 | 1.65% | 3,845.00 | 3,081.29 | 3.37% |
| C9 | 2.93% | 6,644.00 | 5,438.74 | 5.96% |
ENSEIGNEMENT : Le ZONIER présente une fréquence avec beaucoup de variabilité, mais il semble y avoir une croissance par région qu’on peut tester :
| REGION | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| A | 0.92% | 15,169.00 | 13,033.47 | 14.27% |
| B | 1.66% | 32,910.00 | 28,052.53 | 30.72% |
| C | 2.40% | 60,877.00 | 50,227.55 | 55.01% |
3.3.2.10 NBSIN_TYPE1_AN1
| NBSIN_TYPE1_AN1 | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| 0 | 1.83% | 101,793.00 | 84,856.05 | 92.93% |
| 1 | 3.70% | 6,581.00 | 5,945.96 | 6.51% |
| 2 | 4.04% | 530.00 | 469.28 | 0.51% |
| 3 | 9.46% | 52.00 | 42.26 | 0.05% |
ENSEIGNEMENT : La fréquence augmente globalement avec le nombre de sinistres de type 1 l’année précédente. Toutefois, il semble surtout q’il existe un effet d’une indicatrice a eu sinistre ou non.
3.3.2.11 NBSIN_TYPE1_AN2
| NBSIN_TYPE1_AN2 | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| 1 | 1.85% | 102,588.00 | 85,533.45 | 93.67% |
| 2 | 3.50% | 5,822.00 | 5,279.52 | 5.78% |
| 3 | 4.92% | 546.00 | 500.58 | 0.55% |
ATTENTION : comme la table ci-dessous le démontre, cette variable présente un problème puisque toutes les lignes ont au moins un sinistre
Comme discuté dans le hackathon on ignore cette variable
3.3.2.12 NBSIN_TYPE1_AN3
| NBSIN_TYPE1_AN3 | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| 0 | 1.86% | 103,309.00 | 86,105.62 | 94.30% |
| 1 | 3.42% | 5,171.00 | 4,767.33 | 5.22% |
| 2 | 6.13% | 435.00 | 403.13 | 0.44% |
| 3 | 5.34% | 41.00 | 37.46 | 0.04% |
ENSEIGNEMENT : La fréquence augmente globalement avec le nombre de sinistres de type 1 il y a 3 ans. Toutefois, il semble surtout q’il existe un effet d’une indicatrice a eu sinistre ou non.
3.3.2.13 NBSIN_TYPE2_AN1
| NBSIN_TYPE2_AN1 | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| 0 | 1.95% | 106,796.00 | 89,403.29 | 97.91% |
| 1 | 2.61% | 2,034.00 | 1,797.95 | 1.97% |
| 2 | 2.93% | 115.00 | 102.38 | 0.11% |
| 3 | 0.00% | 11.00 | 9.93 | 0.01% |
ENSEIGNEMENT : Il n’existe pas de lien stable sur la fréquence. A ne pas retenir.
3.3.2.14 NBSIN_TYPE2_AN2
| NBSIN_TYPE2_AN2 | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| 0 | 1.96% | 95,382.00 | 79,861.28 | 87.46% |
| 1 | 3.02% | 1,452.00 | 1,307.09 | 1.43% |
| 2 | 1.39% | 77.00 | 71.91 | 0.08% |
| 3 | 0.00% | 3.00 | 3.00 | 0.00% |
| NA | 1.86% | 12,042.00 | 10,070.27 | 11.03% |
ENSEIGNEMENT : Il n’existe pas de lien stable sur la fréquence. A ne pas retenir.
3.3.2.15 NBSIN_TYPE2_AN3
| NBSIN_TYPE2_AN3 | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| 0 | 1.95% | 107,445.00 | 89,931.57 | 98.49% |
| 1 | 2.68% | 1,414.00 | 1,292.12 | 1.42% |
| 2 | 4.31% | 93.00 | 85.86 | 0.09% |
| 3 | 0.00% | 4.00 | 4.00 | 0.00% |
ENSEIGNEMENT : Il n’existe pas de lien stable sur la fréquence. A ne pas retenir.
3.3.2.16 ANNEE
| ANNEE | Freq | nb_lignes | EXPO | pct_EXPO |
|---|---|---|---|---|
| 2016 | 1.82% | 37,947.00 | 31,635.56 | 34.64% |
| 2017 | 2.04% | 36,440.00 | 30,597.37 | 33.51% |
| 2018 | 2.04% | 34,569.00 | 29,080.62 | 31.85% |
ENSEIGNEMENT : L’année 2016 semble atypique. Le niveau actuel semble plus proche de ceux de 2017 et 2018.
3.3.3 Nettoyage des données
preprocess<- function(dt) {
res <- dt %>%
mutate(
Y = NB/EXPO, #calcul de la fréquence pour XGB qui ne supporte pas les offset, mais accepte les poids
NB_PIECES = ifelse(is.na(NB_PIECES), 2, NB_PIECES), #2 est le mode de NB_PIECES
NBSIN_TYPE1_AN1 = ifelse(NBSIN_TYPE1_AN1==0,0,1), #On regroupe les modalités non stables dans le temps
NBSIN_TYPE1_AN3 = ifelse(NBSIN_TYPE1_AN3==0,0,1), #On regroupe les modalités non stables dans le temps
REGION = gsub('[0-9]', '', ZONIER), #création de région
VALEUR_DES_BIENS = pmin(VALEUR_DES_BIENS, 50000), # On regroupe les valeurs des biens supérieurs à 80k
OFFSET = log(EXPO) # création d'une variable d'offset qui servira pour le gbm
) %>%
mutate( REGION = as.factor(REGION),
VALEUR_DES_BIENS = as.factor(VALEUR_DES_BIENS)
) %>%
select(-c('X', 'TYPE_HABITATION', 'NIVEAU_JURIDIQUE', 'NBSIN_TYPE1_AN2', 'NBSIN_TYPE2_AN1', 'NBSIN_TYPE2_AN2', 'NBSIN_TYPE2_AN3', 'ZONIER' )) %>% # On supprime les variables sans lien stable dans le temps
relocate('id', 'EXPO', 'FORMULE', 'TYPE_RESIDENCE', 'SITUATION_JURIDIQUE', 'OBJETS_DE_VALEUR', 'VALEUR_DES_BIENS', 'NB_PIECES', 'REGION', 'NBSIN_TYPE1_AN1', 'NBSIN_TYPE1_AN3', 'ANNEE', 'NB', 'Y', 'COUT' )
return(res)
}
df_train <- preprocess(mrh_train)
df_val <- preprocess(mrh_val)
df_test <- preprocess(mrh_test)3.3.4 Calcul du V de Cramer
On mesude l’association des variables entre elles
X_train <- df_train %>% select(-c('id', 'EXPO', 'NB', 'Y', 'COUT'))
# Initialize empty matrix to store coefficients
empty_m <- matrix(ncol = length(X_train),
nrow = length(X_train),
dimnames = list(names(X_train),
names(X_train)))
# Function that accepts matrix for coefficients and data and returns a correlation matrix
calculate_cramer <- function(m, df) {
for (r in seq(nrow(m))){
for (c in seq(ncol(m))){
m[[r, c]] <- lsr::cramersV(X_train[[r]], X_train[[c]])
}
}
return(m)
}
cor_matrix <- calculate_cramer(empty_m ,X_train)
corrplot(cor_matrix)
4 Modélisation de la fréquence
Dans cette partie, on essaye de modéliser la fréquence des sinistres.
summary(df_test$NB)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01745 0.00000 2.000004.1 Modèle de référence
On commence par une modélisation de référence réduite à un modèle GLM POisson avec un unique intercept.
reg0 = glm(NB~1+offset(OFFSET),family=poisson,data=df_train)
summary(reg0)##
## Call:
## glm(formula = NB ~ 1 + offset(OFFSET), family = poisson, data = df_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.2058 -0.2058 -0.2058 -0.1774 4.8750
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.85523 0.02274 -169.5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 15507 on 108955 degrees of freedom
## Residual deviance: 15507 on 108955 degrees of freedom
## AIC: 19294
##
## Number of Fisher Scoring iterations: 64.2 Modèle GLM plus élaboré
Il s’agit d’un modèle GLM Poisson classique tirant partie des enseignements de la partie précédente, sans interaction. Les variables sont traitées comme catégorielles pour permettre les effets non linéaires.
fit_glm <- glm(NB ~ FORMULE + TYPE_RESIDENCE + SITUATION_JURIDIQUE + OBJETS_DE_VALEUR + as.factor(VALEUR_DES_BIENS) + as.factor(NB_PIECES) + REGION + as.factor(NBSIN_TYPE1_AN1) + as.factor(NBSIN_TYPE1_AN3) + as.factor(ANNEE),
data=df_train, offset=OFFSET, family=quasipoisson())
summary(fit_glm)##
## Call:
## glm(formula = NB ~ FORMULE + TYPE_RESIDENCE + SITUATION_JURIDIQUE +
## OBJETS_DE_VALEUR + as.factor(VALEUR_DES_BIENS) + as.factor(NB_PIECES) +
## REGION + as.factor(NBSIN_TYPE1_AN1) + as.factor(NBSIN_TYPE1_AN3) +
## as.factor(ANNEE), family = quasipoisson(), data = df_train,
## offset = OFFSET)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4898 -0.2195 -0.1769 -0.1220 4.3196
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.35103 0.78732 -10.607 < 2e-16 ***
## FORMULECONFORT 0.89291 0.21816 4.093 4.26e-05 ***
## FORMULEESSENTIEL 0.85630 0.22564 3.795 0.000148 ***
## FORMULEMEDIUM 0.85559 0.21928 3.902 9.55e-05 ***
## TYPE_RESIDENCESECONDAIRE -0.61485 0.11032 -5.573 2.51e-08 ***
## SITUATION_JURIDIQUEPROPRIO -0.05746 0.06009 -0.956 0.338984
## OBJETS_DE_VALEURNIVEAU_2 0.26778 0.08224 3.256 0.001129 **
## as.factor(VALEUR_DES_BIENS)3500 0.62720 0.34944 1.795 0.072676 .
## as.factor(VALEUR_DES_BIENS)9000 0.77173 0.35190 2.193 0.028306 *
## as.factor(VALEUR_DES_BIENS)20000 0.83537 0.35252 2.370 0.017803 *
## as.factor(VALEUR_DES_BIENS)35000 0.97892 0.35803 2.734 0.006254 **
## as.factor(VALEUR_DES_BIENS)50000 0.93078 0.35814 2.599 0.009352 **
## as.factor(NB_PIECES)1 2.04164 0.82753 2.467 0.013621 *
## as.factor(NB_PIECES)2 1.98341 0.82634 2.400 0.016387 *
## as.factor(NB_PIECES)3 2.23041 0.82810 2.693 0.007074 **
## as.factor(NB_PIECES)4 2.31557 0.83137 2.785 0.005350 **
## REGIONB 0.43688 0.10543 4.144 3.42e-05 ***
## REGIONC 0.94976 0.09974 9.522 < 2e-16 ***
## as.factor(NBSIN_TYPE1_AN1)1 0.43828 0.07264 6.034 1.61e-09 ***
## as.factor(NBSIN_TYPE1_AN3)1 0.33205 0.08158 4.070 4.70e-05 ***
## as.factor(ANNEE)2017 0.10332 0.05889 1.754 0.079351 .
## as.factor(ANNEE)2018 0.10273 0.05960 1.724 0.084798 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.119556)
##
## Null deviance: 15507 on 108955 degrees of freedom
## Residual deviance: 14906 on 108934 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 9Le modèle n’explique qu’une faible partie de la variabilité.
Par ailleurs la variable SITUATION_JURIDIQUE (qui n’a que 2 niveaux) n’a finalement pas d’impact sur le risque et peut donc être exclue du modèle :
fit_glm <- glm(NB ~ FORMULE + TYPE_RESIDENCE + OBJETS_DE_VALEUR + as.factor(VALEUR_DES_BIENS) + as.factor(NB_PIECES) + REGION + as.factor(NBSIN_TYPE1_AN1) + as.factor(NBSIN_TYPE1_AN3) + as.factor(ANNEE),
data=df_train, offset=OFFSET, family=quasipoisson())
summary(fit_glm)##
## Call:
## glm(formula = NB ~ FORMULE + TYPE_RESIDENCE + OBJETS_DE_VALEUR +
## as.factor(VALEUR_DES_BIENS) + as.factor(NB_PIECES) + REGION +
## as.factor(NBSIN_TYPE1_AN1) + as.factor(NBSIN_TYPE1_AN3) +
## as.factor(ANNEE), family = quasipoisson(), data = df_train,
## offset = OFFSET)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4903 -0.2202 -0.1776 -0.1224 4.3204
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.41067 0.78465 -10.719 < 2e-16 ***
## FORMULECONFORT 0.90472 0.21774 4.155 3.25e-05 ***
## FORMULEESSENTIEL 0.87116 0.22504 3.871 0.000108 ***
## FORMULEMEDIUM 0.86290 0.21908 3.939 8.20e-05 ***
## TYPE_RESIDENCESECONDAIRE -0.60328 0.10960 -5.504 3.71e-08 ***
## OBJETS_DE_VALEURNIVEAU_2 0.26617 0.08219 3.238 0.001203 **
## as.factor(VALEUR_DES_BIENS)3500 0.61153 0.34890 1.753 0.079655 .
## as.factor(VALEUR_DES_BIENS)9000 0.76481 0.35165 2.175 0.029640 *
## as.factor(VALEUR_DES_BIENS)20000 0.83693 0.35233 2.375 0.017530 *
## as.factor(VALEUR_DES_BIENS)35000 0.98370 0.35781 2.749 0.005974 **
## as.factor(VALEUR_DES_BIENS)50000 0.93646 0.35790 2.616 0.008885 **
## as.factor(NB_PIECES)1 2.06378 0.82697 2.496 0.012577 *
## as.factor(NB_PIECES)2 2.01424 0.82547 2.440 0.014684 *
## as.factor(NB_PIECES)3 2.27101 0.82677 2.747 0.006018 **
## as.factor(NB_PIECES)4 2.35781 0.82997 2.841 0.004500 **
## REGIONB 0.43567 0.10539 4.134 3.57e-05 ***
## REGIONC 0.94502 0.09959 9.489 < 2e-16 ***
## as.factor(NBSIN_TYPE1_AN1)1 0.44183 0.07254 6.091 1.13e-09 ***
## as.factor(NBSIN_TYPE1_AN3)1 0.33622 0.08146 4.128 3.67e-05 ***
## as.factor(ANNEE)2017 0.10379 0.05887 1.763 0.077911 .
## as.factor(ANNEE)2018 0.10351 0.05958 1.737 0.082354 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.119014)
##
## Null deviance: 15507 on 108955 degrees of freedom
## Residual deviance: 14907 on 108935 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 9fit_glm_Diag <- data.frame(df_test,
fit = predict(fit_glm, newdata = df_test, type="response")
) %>%
mutate(pearson = (NB-fit)/sqrt(fit))
dat <- fit_glm_Diag %>%
select( fit, pearson) %>%
mutate(quantile = ntile(fit, 16000)) %>%
group_by(quantile) %>%
summarize(mean_fit = mean(fit, na.rm = TRUE), mean_pearson = mean(pearson, na.rm = TRUE))
ggplot(dat, aes(x = mean_fit, y = mean_pearson)) +
geom_point() +
geom_abline(intercept = 0, slope = 0, size = 1, color="red") +
ggtitle("Rédidus de Pearson aggrégés - test - GLM")
round(sum(fit_glm_Diag$fit)/sum(df_test$EXPO),4)## [1] 0.0214round(sum(df_test$NB)/sum(df_test$EXPO),4)## [1] 0.0209fit_glm_Diag %>%
mutate(residual = NB-fit) %>%
summarize(rmse = sqrt(mean(residual^2)))## rmse
## 1 0.13295754.3 Modèle de Gradient Boosting
Pour la mise en place d’un modèle GBM, nous utiliserons le package h2o :
h2o.init(nthreads = -1)h2o.no_progress()Ce package permet de définir une fonction de perte correspondant à la déviance de Poisson avec sa fonction de lien canonique, c’est à dire le ln dans notre cas. Par ailleurs, la méthode permet de définir une colonne qui peut servir d’offset. D’après la documentation, cet offset est appliqué dans l’espace linéarisé :
For other distributions, the offset corrections are applied in the linearized space before applying the inverse link function to get the actual response values.
Il nous faut donc passer en valeur le \(ln(EXPO)\) que nous avons renseigné dans la variable OFFSET :
En complément nous mettons en oeuvre une recherche par grille pour trouver le métaparamètre de profondeur d’arbre optimal. Le code est adapté de cet article
hyper_params = list( max_depth = seq(1,15,2) )
grid <- h2o.grid(
## hyper parameters
hyper_params = hyper_params,
## full Cartesian hyper-parameter search
search_criteria = list(strategy = "Cartesian"),
## which algorithm to run
algorithm="gbm",
## identifier for the grid, to later retrieve it
grid_id="Grid_gbm_offset",
## standard model parameters
distribution = "poisson",
x = names(df_train)[3:12],
y = 'NB',
offset_column = "OFFSET",
training_frame = as.h2o(df_train),
validation_frame = as.h2o(df_val),
## more trees is better if the learning rate is small enough
## here, use "more than enough" trees - we have early stopping
ntrees = 10000,
## smaller learning rate is better
## since we have learning_rate_annealing, we can afford to start with a bigger learning rate
learn_rate = 0.05,
## learning rate annealing: learning_rate shrinks by 1% after every tree
## (use 1.00 to disable, but then lower the learning_rate)
learn_rate_annealing = 0.99,
## sample 80% of rows per tree
sample_rate = 0.8,
## sample 80% of columns per split
col_sample_rate = 0.8,
## fix a random number generator seed for reproducibility
seed = 1234,
## early stopping once the validation AUC doesn't improve by at least 0.01% for 5 consecutive scoring events
stopping_rounds = 5,
stopping_tolerance = 1e-4,
stopping_metric = "deviance",
## score every 10 trees to make early stopping reproducible (it depends on the scoring interval)
score_tree_interval = 10
)
## by default, display the grid search results sorted by increasing logloss (since this is a classification task)
grid## H2O Grid Details
## ================
##
## Grid ID: Grid_gbm_offset
## Used hyper parameters:
## - max_depth
## Number of models: 8
## Number of failed models: 0
##
## Hyper-Parameter Search Summary: ordered by increasing residual_deviance
## max_depth model_ids residual_deviance
## 1 3 Grid_gbm_offset_model_2 0.17448113788335162
## 2 1 Grid_gbm_offset_model_1 0.17561652746681825
## 3 5 Grid_gbm_offset_model_3 0.17986113859022135
## 4 7 Grid_gbm_offset_model_4 0.19916164139992715
## 5 9 Grid_gbm_offset_model_5 0.24685830723796695
## 6 11 Grid_gbm_offset_model_6 0.3093800261468491
## 7 13 Grid_gbm_offset_model_7 0.3404180725614419
## 8 15 Grid_gbm_offset_model_8 0.34838731404420437On récupère le meilleur modèle :
fit_gbm <- h2o.getModel(grid@model_ids[[1]])On vérifie la selection en comparant les performances sur le jeu de validation par rapport aux valeurs affichées dans la grille :
print(h2o.performance(fit_gbm, newdata = as.h2o(df_val)))## H2ORegressionMetrics: gbm
##
## MSE: 0.01898311
## RMSE: 0.1377792
## MAE: 0.03495651
## RMSLE: 0.0935029
## Mean Residual Deviance : 0.1744811plot(fit_gbm)
h2o.performance(fit_gbm, newdata = as.h2o(df_test))## H2ORegressionMetrics: gbm
##
## MSE: 0.0176801
## RMSE: 0.1329665
## MAE: 0.03420818
## RMSLE: 0.0912325
## Mean Residual Deviance : 0.1701883fit_gbm_Diag <- data.frame(df_test,
fit = as.vector(h2o.predict(object = fit_gbm, newdata=as.h2o(df_test)))
) %>%
mutate(pearson = (NB-fit)/sqrt(fit))
dat_gbm <- fit_gbm_Diag %>%
select( fit, pearson) %>%
mutate(quantile = ntile(fit, 16000)) %>%
group_by(quantile) %>%
summarize(mean_fit = mean(fit, na.rm = TRUE), mean_pearson = mean(pearson, na.rm = TRUE))
ggplot(dat_gbm, aes(x = mean_fit, y = mean_pearson)) +
geom_point() +
geom_abline(intercept = 0, slope = 0, size = 1, color="red") +
ggtitle("Rédidus de Pearson aggrégés - test - GBM")
fit_gbm_Diag %>%
mutate(residual = NB-fit) %>%
summarize(rmse = sqrt(mean(residual^2)))## rmse
## 1 0.13296654.4 Modèle de XGBoost
x <- c('FORMULE', 'TYPE_RESIDENCE', 'SITUATION_JURIDIQUE','OBJETS_DE_VALEUR', 'VALEUR_DES_BIENS', 'NB_PIECES', 'REGION', 'NBSIN_TYPE1_AN1', 'NBSIN_TYPE1_AN3', 'ANNEE')
y <- 'Y'
w <- 'EXPO'# Input maker
prep_xgb <- function(dat, x) {
data.matrix(dat[, x, drop = FALSE])
}
# Data interface to XGBoost
dtrain <- xgb.DMatrix(
prep_xgb(df_train, x),
label = df_train[[y]],
weight = df_train[[w]]
)
# Parameters chosen by 5-fold grouped CV
params_freq <- list(
learning_rate = 0.2,
max_depth = 5,
alpha = 3,
lambda = 0.5,
max_delta_step = 2,
min_split_loss = 0,
colsample_bytree = 1,
subsample = 0.9
)
# Fit
set.seed(1)
fit_xgb <- xgb.train(
params_freq,
data = dtrain,
nrounds = 580,
objective = "count:poisson",
watchlist = list(train = dtrain),
print_every_n = 100
)## [1] train-poisson-nloglik:0.504203
## [101] train-poisson-nloglik:0.114903
## [201] train-poisson-nloglik:0.100885
## [301] train-poisson-nloglik:0.100260
## [401] train-poisson-nloglik:0.099847
## [501] train-poisson-nloglik:0.099478
## [580] train-poisson-nloglik:0.099240# Save and load model
xgb.save(fit_xgb, "xgb.model")## [1] TRUEfit_xgb <- xgb.load("xgb.model")fit_xgb_Diag <- data.frame(df_test,
fit = predict(fit_xgb, prep_xgb(df_test, x))
) %>%
mutate(pearson = (NB-fit)/sqrt(fit))
dat_xgb <- fit_xgb_Diag %>%
select( fit, pearson) %>%
mutate(quantile = ntile(fit, 16000)) %>%
group_by(quantile) %>%
summarize(mean_fit = mean(fit, na.rm = TRUE), mean_pearson = mean(pearson, na.rm = TRUE))
ggplot(dat_xgb, aes(x = mean_fit, y = mean_pearson)) +
geom_point() +
geom_abline(intercept = 0, slope = 0, size = 1, color="red") +
ggtitle("Rédidus de Pearson aggrégés - test - XGB")
fit_xgb_Diag %>%
mutate(residual = NB-fit) %>%
summarize(rmse = sqrt(mean(residual^2)))## rmse
## 1 0.13309365 Explicabilité des modèles produits
Cette section se consacre à l’explicabilité des modèles produits dans le but d’identifier la manière dont les modèles ML exploitent les combinaisons de variables et effets non linéaires.
Cette section s’appuie sur les travaux présentés dans Peeking into the black box de Christian Lorentzen et Michael Mayer
5.1 Modèles de prédiction
La première étape consiste à définir des objets génériques pour obtenir des prédiction des modèles. Nous allons nous appuyer sur le package flashlight en définissant un objet par modèle, mais qui aura par la suite le même structure d’appel.
fl_glm <- flashlight(
model = fit_glm, label = "GLM",
predict_function = function(fit, X) predict(fit, X, type = "response")
)
fl_xgb <- flashlight(
model = fit_xgb, label = "XGBoost",
predict_function = function(fit, X) predict(fit, prep_xgb(X, x))
)
fl_gbm <- flashlight(
model = fit_gbm, label = "H20 GBM",
predict_function = function(fit, X) as.vector(h2o.predict(object = fit_gbm, newdata=as.h2o(X)))
)
# On combine les différents modèles dans un objet unique disposant des métriques définies ci-dessous
metrics <- list(
`Déviance Poisson` = deviance_poisson,
`R² Déviance Poisson` = r_squared_poisson)
fls <- multiflashlight(list(fl_glm, fl_xgb, fl_gbm), data = df_test,
y = y, w = w, metrics = metrics)
# Variation de l'objet précédent dans l'espace linéarisé
fls_log <- multiflashlight(fls, linkinv = log)5.2 Comparaison de la performance des modèles
On commence par comparer la performance globale des modèles ajustés jusqu’à maintenant :
fillc <- "#E69F00"
perf <- light_performance(fls)
perf##
## I am an object of class light_performance_multi
##
## data.frames:
##
## data
## # A tibble: 6 x 3
## label metric value
## <fct> <fct> <dbl>
## 1 GLM Déviance Poisson 0.166
## 2 GLM R² Déviance Poisson 0.0107
## 3 XGBoost Déviance Poisson 0.163
## 4 XGBoost R² Déviance Poisson 0.0280
## 5 H20 GBM Déviance Poisson 0.166
## 6 H20 GBM R² Déviance Poisson 0.0103plot(perf, geom = "point") +
labs(x = element_blank(), y = element_blank())
Du point de vue de la déviance de Poisson, les différents modèles présentent des performances comparables avec un léger avantage au modèke XGBoost. C’est un des éléments qui était préssenti : le jeu de donné ne semble pas permettre une mise en avant des méthodes non linéaires. Un des poinst saillants est la faible fréquence moyenne observée de l’orde de 2,12% soit un écart type de 14,55%.
sprintf('férquence moyenne : %.2f%%', tx*100)## [1] "férquence moyenne : 2.12%"sprintf('férquence écart type : %.2f%%', sqrt(tx)*100)## [1] "férquence écart type : 14.55%"Cela explique en partie la faible part de variabilité expliquée par les modèles.
5.3 Importance
On peut ensuite regarder l’importance (au sens permutation) des variables dans chaque modèle. Une version de décomposition de Shapley est disponible mais trop lente pour les moyens à disposition pour ce projet.
imp <- light_importance(fls, v = x, type = 'permutation')
plot(imp, fill = fillc, color = "black")
5.3.1 Partial dependence curves
Pour comprendre la manière dont les différentes varoables contribuent en moyennes au prédiction des modèles on peut s’appuyer sur les pdp. On représente ici les pdp des 3 variables les plus importantes à priori issues de la section précédente :
plot(light_profile(fls, v = "REGION"))
plot(light_profile(fls, v = "VALEUR_DES_BIENS"))
plot(light_profile(fls, v = "FORMULE"))
On peut remarquer que si de légères différences existent dans le niveau relatif de chaque modalité, les 3 modèles sont globalement en accord sur la forme des réponses.
5.3.2 Comparaison des différences de prise en compte de l’effet de la variable VALEUR_DES_BIENS
Dans le graphique ci-dessous on superpose le tri à plat des réponses observées, des prédictions du modèle (en prenant toutes les variables en considération) et les pdp. Ce graphique donne une intuition de la manière dont un variable contribue en moyenne aux performance du modèle et en fonction des zones sur lesquelles il existe plus ou moins d’exposition. On y observe par exemple que la réponse combinée du XGBoost est la plus proche de la réponse observée (comme attendu par rapport aux performances globales du modèle).
eff_VALEUR_DES_BIENS <- light_effects(fls, v = "VALEUR_DES_BIENS", counts_weighted = TRUE)
p <- plot(eff_VALEUR_DES_BIENS, show_points = FALSE)
plot_counts(p, eff_VALEUR_DES_BIENS, alpha = 0.3)
Ce qui est plus surprenant c’est le désaccord qu’ont les modèles sur l’effet des biens de 3500 Eur, qui est portant la modalité qui dispose du plus d’exposition. Le XGBoost y voit :
un effet important plus importants que les autres modèles
un effet relatif d’un bien à 3500 Eur bien plus important que celui d’un bien de valeur nulle.
5.3.3 Analyse des intéractions
Un des points délicats de l’analyse des GLM est l’identification des intéractions potentielles à prendre en compte, là où les méthodes à base d’arbre n’ont pas cette difficulté. Le graphique suivant se base se le H de Friedman qui calcule l’effet de la pdp croisée nette des effets multipliés des pdp univariées :
interact_rel <- light_interaction(
fls_log,
v = most_important(imp, 4),
take_sqrt = FALSE,
pairwise = TRUE,
use_linkinv = TRUE,
seed = 61
)
plot(interact_rel, color = "black", fill = fillc, rotate_x = TRUE)
On y observe que (par construction) le modèle GLM n’inclut aucune intéraction, là où les modèles XGBoost et GBM en considèrent. Toutefois les deux modèles d’arbres n’ont pas la même vision de l’importance de ces intéractions. Si les deux s’accordent sur l’existence d’une intéraction entre VALEUR_DES_BIENS et TYPE_RESIDENCE, XGBoost accorde par ailleurs de l’importance à l’intéraction VALEUR_DES_BIENS et NB_PIECES, là où GBM y voit un intérêt marginal. Les effets ici sont normalisés. On peut aussi les représenter dans leur valeur absolue :
interact_abs <- light_interaction(
fls_log,
v = most_important(imp, 4),
normalize = FALSE,
pairwise = TRUE,
use_linkinv = TRUE,
seed = 61
)
plot(interact_abs, color = "black", fill = fillc, rotate_x = TRUE)
5.3.4 Visualisation
Pour essayer de comprendre la forme de ces intéractions, on peut représenter les pdp d’une variable en fonction des modalités d’une autre :
#Intéractions fortes
pdp_VALEUR_DES_BIENS_NB_PIECES <- light_profile(fls_log, v = "VALEUR_DES_BIENS", by = "NB_PIECES",
pd_seed = 50, data = df_val, counts_weighted=TRUE)
plot(pdp_VALEUR_DES_BIENS_NB_PIECES)
pdp_VALEURDESBIENS_TYPE_RESIDENCE <- light_profile(fls_log, v = "VALEUR_DES_BIENS", by = "TYPE_RESIDENCE",
pd_seed = 50, data = df_val, counts_weighted=TRUE)
plot(pdp_VALEURDESBIENS_TYPE_RESIDENCE)
# Intéractions faibles
pdp_TYPE_RESIDENCE_NB_PIECES <- light_profile(fls_log, v = "TYPE_RESIDENCE",
by = "NB_PIECES", pd_seed = 50, data = df_val, counts_weighted=TRUE)
plot(pdp_TYPE_RESIDENCE_NB_PIECES)
Le premier graphique entre VALEUR_DES_BIENS et NB_PIECES est très intéressant et illustre à la fois une divergence de vue entre les modèles et explique la forme de la réponse moyenne de la VALER_DES_BIENS du XGBoost précédement vue.
En effet tous les modèles s’accordent sur une croissance du risque avec le nombre de pièces. Par ailleurs, le GLM comme le XGBoost s’accordent sur un risque très inférieur sur les biens de valeur nulle sans pièces (c’est la limite de cet exercice sur des données anonymisées, il n’est pas possible de revenir auprès du métier pour comprendre le rationnel de ce codage). Par contre, le GLM, sans intéraction, se voit contraint d’accorder un risque homogène en fonction de la valeur du bien, là où le XGBoost, plus libre accorde un niveau de risque beaucoup plsu important aux biens sans pièce, mais de valeur non nulle. En outre le XGBoost comme le GBM dans une moindre mesure sur une décroissance du risque des biens de valeur intermédiaire.
En complément XGBoost comme GBM s’accordent sur un plafonnement du risque pour les résidences secondaires en fonction de la valeur des biens (contrairement aux résidences principales).
6 Apport dans la modélisation GLM
Fort des constats précédents on peut créer un nouveau GLM incluant une intéraction entre la VALEUR_DES_BIENS et le NB_PIECES. Cette étape nécessite des allers retours en fonction de l’évolution de la pertinence des variabes dans le modèle suite à cette adjunction. Idéalement il aurait fallu créer une nouvelle variable NB_PIECES avec 2 modalités 0 et 1+ et créer une intéraction avec cette dernière pour recréer l’effet plateau sans augmenter de manière indue le nombre de paramètres du modèle.
fit_glm_opti <- glm(NB ~ FORMULE + TYPE_RESIDENCE + OBJETS_DE_VALEUR + as.factor(VALEUR_DES_BIENS)+ REGION + as.factor(NBSIN_TYPE1_AN1) + as.factor(NBSIN_TYPE1_AN3) + as.factor(ANNEE) + as.factor(VALEUR_DES_BIENS):NB_PIECES,
data=df_train, offset=OFFSET, family=quasipoisson())
summary(fit_glm_opti)##
## Call:
## glm(formula = NB ~ FORMULE + TYPE_RESIDENCE + OBJETS_DE_VALEUR +
## as.factor(VALEUR_DES_BIENS) + REGION + as.factor(NBSIN_TYPE1_AN1) +
## as.factor(NBSIN_TYPE1_AN3) + as.factor(ANNEE) + as.factor(VALEUR_DES_BIENS):NB_PIECES,
## family = quasipoisson(), data = df_train, offset = OFFSET)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5010 -0.2203 -0.1805 -0.1237 4.3424
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.60321 0.67661 -12.715 < 2e-16
## FORMULECONFORT 0.90899 0.21780 4.173 3.00e-05
## FORMULEESSENTIEL 0.87623 0.22509 3.893 9.92e-05
## FORMULEMEDIUM 0.86441 0.21921 3.943 8.04e-05
## TYPE_RESIDENCESECONDAIRE -0.58754 0.10945 -5.368 7.98e-08
## OBJETS_DE_VALEURNIVEAU_2 0.26529 0.08214 3.230 0.001240
## as.factor(VALEUR_DES_BIENS)3500 2.72191 0.64697 4.207 2.59e-05
## as.factor(VALEUR_DES_BIENS)9000 2.77466 0.65256 4.252 2.12e-05
## as.factor(VALEUR_DES_BIENS)20000 3.12283 0.65341 4.779 1.76e-06
## as.factor(VALEUR_DES_BIENS)35000 2.67128 0.68710 3.888 0.000101
## as.factor(VALEUR_DES_BIENS)50000 2.81803 0.67656 4.165 3.11e-05
## REGIONB 0.44238 0.10543 4.196 2.72e-05
## REGIONC 0.94344 0.09958 9.474 < 2e-16
## as.factor(NBSIN_TYPE1_AN1)1 0.44137 0.07254 6.085 1.17e-09
## as.factor(NBSIN_TYPE1_AN3)1 0.33627 0.08146 4.128 3.66e-05
## as.factor(ANNEE)2017 0.10322 0.05887 1.753 0.079539
## as.factor(ANNEE)2018 0.10402 0.05958 1.746 0.080809
## as.factor(VALEUR_DES_BIENS)0:NB_PIECES 1.03680 0.24749 4.189 2.80e-05
## as.factor(VALEUR_DES_BIENS)3500:NB_PIECES 0.08672 0.07857 1.104 0.269753
## as.factor(VALEUR_DES_BIENS)9000:NB_PIECES 0.12851 0.07367 1.744 0.081103
## as.factor(VALEUR_DES_BIENS)20000:NB_PIECES 0.01112 0.06517 0.171 0.864533
## as.factor(VALEUR_DES_BIENS)35000:NB_PIECES 0.24970 0.09164 2.725 0.006435
## as.factor(VALEUR_DES_BIENS)50000:NB_PIECES 0.17756 0.07515 2.363 0.018147
##
## (Intercept) ***
## FORMULECONFORT ***
## FORMULEESSENTIEL ***
## FORMULEMEDIUM ***
## TYPE_RESIDENCESECONDAIRE ***
## OBJETS_DE_VALEURNIVEAU_2 **
## as.factor(VALEUR_DES_BIENS)3500 ***
## as.factor(VALEUR_DES_BIENS)9000 ***
## as.factor(VALEUR_DES_BIENS)20000 ***
## as.factor(VALEUR_DES_BIENS)35000 ***
## as.factor(VALEUR_DES_BIENS)50000 ***
## REGIONB ***
## REGIONC ***
## as.factor(NBSIN_TYPE1_AN1)1 ***
## as.factor(NBSIN_TYPE1_AN3)1 ***
## as.factor(ANNEE)2017 .
## as.factor(ANNEE)2018 .
## as.factor(VALEUR_DES_BIENS)0:NB_PIECES ***
## as.factor(VALEUR_DES_BIENS)3500:NB_PIECES
## as.factor(VALEUR_DES_BIENS)9000:NB_PIECES .
## as.factor(VALEUR_DES_BIENS)20000:NB_PIECES
## as.factor(VALEUR_DES_BIENS)35000:NB_PIECES **
## as.factor(VALEUR_DES_BIENS)50000:NB_PIECES *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.118798)
##
## Null deviance: 15507 on 108955 degrees of freedom
## Residual deviance: 14903 on 108933 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 9On peut alors rajouter ce modèle à notre corpus et en comparer la performance :
fl_glm_opti <- flashlight(
model = fit_glm_opti, label = "GLM_opti",
predict_function = function(fit, X) predict(fit, X, type = "response")
)
# On combine les différents modèles dans un objet unique disposant des métriques définies ci-dessous
fls_opti <- multiflashlight(list(fl_glm, fl_glm_opti, fl_xgb, fl_gbm), data = df_test,
y = y, w = w, metrics = metrics)
# Variation de l'objet précédent dans l'espace linéarisé
fls_log_opti <- multiflashlight(fls, linkinv = log)fillc <- "#E69F00"
perf_opti <- light_performance(fls_opti)
perf_opti##
## I am an object of class light_performance_multi
##
## data.frames:
##
## data
## # A tibble: 8 x 3
## label metric value
## <fct> <fct> <dbl>
## 1 GLM Déviance Poisson 0.166
## 2 GLM R² Déviance Poisson 0.0107
## 3 GLM_opti Déviance Poisson 0.166
## 4 GLM_opti R² Déviance Poisson 0.0118
## 5 XGBoost Déviance Poisson 0.163
## 6 XGBoost R² Déviance Poisson 0.0280
## 7 H20 GBM Déviance Poisson 0.166
## 8 H20 GBM R² Déviance Poisson 0.0103plot(perf_opti, geom = "point") +
labs(x = element_blank(), y = element_blank())
Le gain ici reste marginal et mériterait d’être encore travaillé par rapport aux autres intéractions détectées dans les modèles, néanmoins elle illustre la manière dont les enseignements des modèles ML peuvent être explicités et servir de guide dans l’amélioration de modèles GLM.
7 Conclusion
Ce TD a été l’occasion de poursuivre le travail entrepris lors du Hackathon. Ce jeu de données n’était peut-être pas le plus adapté du fait de la simplicité de sa structure et de la faible fréquence et donc l’importance de la variabilité du phénomène suivi. Néanmoins il nous a permi de remprendre la démarche présentée dans Peeking into the black box pour illustrer comment les approches de ML peuvent utilement aiguiller le praticien dans la construction des GLM qui sont prédominant dans les études tarifaires en entreprise.